cleanup before you start
rm( list=ls() )
library(moments)
setup
set.seed(1234)
N <- 500
H <- 5
phi0 <- 0.2
phi1 <- 0.97
theta1 <- 0.6
sig.eps <- 0.1
simulate DGP
y <- rep(0,N+H)
eps <- rnorm(N+H,0,sig.eps)
y[1] <- phi0/(1.0-phi1) + eps[1]
for( t in 2:(N+H) ){
y[t] <- phi0 + phi1*y[t-1] + eps[t] - theta1*eps[t-1]
}
y.out <- y[(N+1):(N+H)]
y <- y[1:N]
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
grid( lwd=1 , col="darkgrey" )
LBQ Test
Box.test( y, lag=22 , type="Ljung-Box" )
##
## Box-Ljung test
##
## data: y
## X-squared = 4597.3, df = 22, p-value < 2.2e-16
# ACF & PACF
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( y , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
Estimate ARMA Models
ar1 = arima(y,order=c(1,0,0))
ar2 = arima(y,order=c(2,0,0))
ma1 = arima(y,order=c(0,0,1))
arma11 = arima(y,order=c(1,0,1))
ar1_aic <- (-2*ar1$loglik+2*3)/N
ar2_aic <- (-2*ar2$loglik+2*4)/N
ma1_aic <- (-2*ma1$loglik+2*3)/N
arma11_aic <- (-2*arma11$loglik+2*4)/N
ar1_bic <- (-2*ar1$loglik+log(N)*3)/N
ar2_bic <- (-2*ar2$loglik+log(N)*4)/N
ma1_bic <- (-2*ma1$loglik+log(N)*3)/N
arma11_bic <- (-2*arma11$loglik+log(N)*4)/N
round( rbind( c(ar1$loglik,ar2$loglik,ma1$loglik,arma11$loglik),
c(ar1_aic,ar2_aic,ma1_aic,arma11_aic) ,
c(ar1_bic,ar2_bic,ma1_bic,arma11_bic) ) , 3 )
## [,1] [,2] [,3] [,4]
## [1,] 383.158 410.141 234.417 425.704
## [2,] -1.521 -1.625 -0.926 -1.687
## [3,] -1.495 -1.591 -0.900 -1.653
FITTED VALUES
ar1_mu <- y-ar1$residuals
ar2_mu <- y-ar2$residuals
ma1_mu <- y-ma1$residuals
arma11_mu <- y-arma11$residuals
ar1_res <- as.numeric(ar1$residuals)
ar2_res <- as.numeric(ar2$residuals)
ma1_res <- as.numeric(ma1$residuals)
arma11_res <- as.numeric(arma11$residuals)
Plot a bunch of stuff..
# AR1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ar1_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ar1_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ar1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# AR2
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ar2_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ar2_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ar2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ar2_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# MA1
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( ma1_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( ma1_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( ma1_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
# ARMA11
par( mar=c(2,2,1,1) , xaxs="i" , mfrow=c(2,1) )
plot( y , t='p' , pch=16, col='darkorange2' , tck = 0.02 , xlim=c(1,N+H) )
lines( arma11_mu , t='l' , lwd=2 , col='blue3' )
grid( lwd=1 , col="darkgrey" )
plot( arma11_res , col='purple' )
abline( h=0 , lwd=2 )
grid( lwd=1 , col="darkgrey" )
par( mar=c(2,2,1,1) , mfrow=c(2,1) )
acf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('ACF'),col=c('darkorange2'),lwd=3)
pacf( arma11_res , ylim=c(-0.2,1) , lwd=5 , xlim=c(0,25) , col='darkorange2' , tck=0.02)
legend('topright',c('PACF'),col=c('darkorange2'),lwd=3)
Residuals
par( mar=c(2,2,3,2) , mfrow=c(2,2) )
kernel <- density(ar1_res/sqrt(ar1$sigma2))
plot( kernel , main='AR1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(ar2_res/sqrt(ar2$sigma2))
plot( kernel , main='AR2' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(ma1_res/sqrt(ma1$sigma2))
plot( kernel , main='MA1' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
kernel <- density(arma11_res/sqrt(arma11$sigma2))
plot( kernel , main='ARMA11' )
polygon( kernel , col="tomato" , border='darkred')
abline(h=0,lwd=2)
lines( seq(-10,20,0.1) , dnorm( seq(-10,20,0.1) ) , col='darkblue' ,lwd=2 )
par( mar=c(2,2,3,2) , mfrow=c(2,2) )
qqnorm(ar1_res,col='tomato',main='AR1')
qqline(ar1_res,lwd=2,lty=3)
qqnorm(ar2_res,col='tomato',main='AR2')
qqline(ar2_res,lwd=2,lty=3)
qqnorm(ma1_res,col='tomato',main='MA1')
qqline(ma1_res,lwd=2,lty=3)
qqnorm(arma11_res,col='tomato',main='ARMA11')
qqline(arma11_res,lwd=2,lty=3)
Residual Test
ar1.res.jb <- jarque.test(ar1_res)$p.value
ar2.res.jb <- jarque.test(ar2_res)$p.value
ma1.res.jb <- jarque.test(ma1_res)$p.value
arma11.res.jb <- jarque.test(arma11_res)$p.value
ar1.res.lbq <- Box.test( ar1_res, lag=22 , type="Ljung-Box" )$p.value
ar2.res.lbq <- Box.test( ar2_res, lag=22 , type="Ljung-Box" )$p.value
ma1.res.lbq <- Box.test( ma1_res, lag=22 , type="Ljung-Box" )$p.value
arma11.res.lbq <- Box.test( arma11_res, lag=22 , type="Ljung-Box" )$p.value
round( rbind( c(ar1.res.jb,ar2.res.jb,ma1.res.jb,arma11.res.jb) ,
c(ar1.res.lbq,ar2.res.lbq,ma1.res.lbq,arma11.res.lbq) ) , 3 )
## [,1] [,2] [,3] [,4]
## [1,] 0.123 0.436 0.001 0.424
## [2,] 0.000 0.000 0.000 0.452
FORECAST
ar1_pred <- predict( ar1 , n.ahead=H )
ar2_pred <- predict( ar2 , n.ahead=H )
ma1_pred <- predict( ma1 , n.ahead=H )
arma11_pred <- predict( arma11 , n.ahead=H )
ar1_mse = mean( (y.out - as.numeric(ar1_pred$pred) )**2 )
ar2_mse = mean( (y.out - as.numeric(ar2_pred$pred) )**2 )
ma1_mse = mean( (y.out - as.numeric(ma1_pred$pred) )**2 )
arma11_mse = mean( (y.out - as.numeric(arma11_pred$pred) )**2 )
plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('AR(1) MSE %3.3f',ar1_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange' )
abline( v=N , lwd=2 )
abline( h=ar1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ar1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ar1_pred$pred) , t='b' , lwd=2 , col='blue3' )
plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('AR(2) MSE %3.3f',ar2_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange')
abline( v=N , lwd=2 )
abline( h=ar2$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ar2_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ar2_pred$pred) , t='b' , lwd=2 , col='blue3' )
plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('MA(1) MSE %3.3f',ma1_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange')
abline( v=N , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , ma1_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(ma1_pred$pred) , t='b' , lwd=2 , col='blue3' )
plot( c((N-10):(N+H)) , c(y[(N-10):N], y.out) , main=sprintf('ARMA(1,1) MSE %3.3f',arma11_mse) , ylim=c(min(y),max(y)) , ylab='',xlab='', tck = 0.02 , pch=16 , col='darkorange')
abline( v=N , lwd=2 )
abline( h=ma1$coef['intercept'] , lwd=2 )
grid( lwd=1 , col="darkgrey" )
lines( c((N-10):N) , arma11_mu[(N-10):N] , t='l' , lwd=2 , col='blue3' )
lines( c((N+1):(N+H)) , as.numeric(arma11_pred$pred) , t='b' , lwd=2 , col='blue3' )
round(c( ar1_mse , ar2_mse , ma1_mse , arma11_mse )*100 , 3 )
## [1] 1.381 1.104 4.084 1.043